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Abstract 

The analysis of electromagnetic scattering has long been performed on a discrete representation 
of the geometry. This representation is typically continuous but not differentiable. The need to 
define physical quantities on this geometric representation has led to development of sets of ba¬ 
sis functions that need to satisfy constraints at the boundaries of the elements/tesselations (viz., 
continuity of normal or tangential components across element boundaries). For electromagnetics, 
these result in either curl/div-conforming basis sets. The geometric representation used for analysis 
is in stark contrast with that used for design, wherein the surface representation is higher order 
differentiable. Using this representation for both geometry and physics on geometry has several ad¬ 
vantages, and is eludicated in Hughes et ah. Isogeometric analysis: CAD, finite elements, NURBS, 
exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering 
194 (39-41) (2005). Until now, a bulk of the literature on isogeometric methods have been limited 
to solid mechanics, with some effort to create NURBS based basis functions for electromagnetic 
analysis. In this paper, we present the first complete isogeometry solution methodology for the 
electric field integral equation as applied to simply connected structures. This paper systematically 
proceeds through surface representation using subdivision, definition of vector basis functions on 
this surface, to fidelity in the solution of integral equations. We also present techniques to stabilize 
the solution at low frequencies, and impose a Calderon preconditioner. Several results presented 
serve to validate the proposed approach as well as demonstrate some of its capabilities. 

Keywords: Isogeometric Analysis, Electric Eield Integral Equation, Calderon Identity, 
Low-Erequency Breakdown, Electromagnetics 


1. Introduction 

Computational methods have become the mainstay of scientific investigation in numerous dis¬ 
ciplines, and electromagnetics is no exception. Research in both integral equation and differential 
equation based methods has grown by leaps and bounds over the past few decades. This period 
has witnessed development of both higher order basis functions [H El El n E], and higher order 
representations of geometry (at least locally) O IH |6] amongst many other equally important ad¬ 
vances. These approaches have been applied to a wide range of realistic problems spanning several 
wavelengths. 

However, despite advances in these areas, there is a fundamental disconnect between the ge¬ 
ometry processing and analysis based on this geometry. Traditional analysis proceeds by defining 
a discrete representation of the geometry typically comprising piecewise continuous tessellations. 
Ironically, this discrete representation of the geometry is obtained using software or a computer 
aided design (CAD) tool that contains a higher order differentiable representation of the geometry. 
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As eloquently elucidated in [7], the rationale for this disconnect can be attributed to the different 
periods in time that CAD tools and analysis tools developed. As the latter is older, the computa¬ 
tional foundation is older as well. As a result, one is left with awkward communication with the 
CAD software for refining and remeshing. This is especially true insofar as accuracy is concerned; 
lack of higher order continuity in geometry can cause artifacts if the underlying spaces for field 
representations are not properly defined. Indeed, the need to define div/curl conforming spaces 
on tessellations that are only Co led to development of novel basis sets that meet this criterion 
[8]. An alternate approach that has recently been espoused is isogeometric analysis (IGA). In 
this approach, the basis functions used to represent the geometry are the same as those used to 
represent the physics on this geometry. As a result, the features of geometry representation such 
as higher order continuity, adaptivity, etc., carry over to function representation as well. Since 
the appearance of this approach in 2005 [7], it has been applied to a number of applications that 
range from structure mechanics [9] to fluid-structure interactions (FSI) [1^ to contact problems 
m to flow [12] to shell analysis uadi] to acoustics [15] and electromagnetics m- In addition 
to analysis techniques, the power of IGA has been harnessed for design-through-analysis phase in 
several practical applications [I71IIH1[IS]. 

Next, we briefly review some of the existing methods. Most CAD tools use bi/tri-variate 
spline based patches/solids like those based upon Bezier, B-splines, and non-uniform rational B- 
splines (NURBS). As a result, these basis functions are the most often used as IGA basis, with 
the most popular being NURBS. The latter choice is determined by the fact that NURBS is the 
industry standard for modern CAD systems. Properties such as non-negativity and the fact that 
it provides a partition of unity make it an excellent candidate for defining function spaces. Finite 
element methods based on NURBS basis functions that exhibit h- and p-adaptivity have been 
demonstrated [7]. Unfortunately, the challenge with using NURBS arises from the fact that the 
resulting shapes are topologically either a disk, a tube or a torus. As a result, stitching together 
these patches can result in surfaces that are not watertight. These complexities are exacerbated 
when the object being meshed is topologically complex or has multiple scales [101 HI]- Two other 
geometry processing methodologies gaining popularity for handling shapes that are complex are T- 
splines and subdivision surface. The former, an extension to NURBS, can handle T-junctions and 
hence and greatly reduce the number of the control points in the control mesh. T-splines, especially 
analysis-ready T-splines, comprise a good candidate for constructing isogeometric analysis. More 
detailed work on T-splines and its application in IGA can be found in [n m] and references 
therein. As opposed to T-splines, subdivision surfaces have played a significant role in the computer 
animation industry. Among its many advantages are the ease with which one can represent complex 
topologies, scalability, inherently multiresolution features, efficiency and ease of implementation. 
Furthermore, it converges to a smooth limit surface that is almost everywhere except at isolated 
points where it is continuous. There are several subdivision schemes; Loop, Gatmull-Glark, Doo- 
Sabin to name a few. Generally speaking, all of the three of these schemes alluded above can be 
used to construct IGA method. To date, isogeometric analysis based on subdivision surface is less 
well studied. Some work on IGA based on Catmull-Clark can be found in [23], where IGA is used 
to solve PDFs defined on a surface. 

While the literature on IGA for differential equations is reasonably widespread across multiple 
fields, IGA for integral equations is still at a nascent stage. As a result, it has recently become the 
focus of significant attention. Recently, two dimensional isogeometric boundary element method 
(IGBEM) was proposed [24] to study elastostatic problem with NURBS interpolating basis to rep¬ 
resent the geometry, displacement and tractions. In [25], IGBEM based on unstructured T-splines 
was developed for a three dimensional linear elastostatic problem. This approach was extended to 
address lEs associated with hydrodynamic interactions m- Likewise, IGBEM methods have been 
developed to study acoustic scattering from rigid bodies [IS] as well as two-dimensional electro¬ 
magnetic analysis [2^. To our knowledge, IGA has not been used for solution to three dimensional 
lEs associated with vector electromagnetic fields, and this serves to motivate this paper. 

The focus of this paper will be on the construction of a well formulated low-frequency stable 
IGA solver for the electric field integral equation (EFIE) that is based on subdivision surfaces 
specifically, the Loop subdivision scheme. As will be evident, the choice of Loop subdivision scheme 
is only incidental; the presented method can be applied to most subdivision surface description. 
In developing a solver that is robust, several challenges need to be addressed; these range from 
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definition of basis functions that correctly map the trace of fields on the surface to formulations 
that render the resulting system frequency stable to formulation of effective preconditioners. To 
set the stage for introduction of this formulation, we will assume that the surfaces are simply 
connected and have smoothness almost everywhere (extension to surfaces with sharp edges and 
corners is under development and will form the basis of a subsequent paper). As will be evident, 
the assumption of sufficient “local” smoothness permits significant freedom in terms of defining 
function spaces. Thus, the principal contribution of this work is four-fold: We will 

• present construction of a basis to correctly represent the trace fields on surface with genus 
equal to zero, 

• demonstrate convergence of EFIE-IGA solver for canonical geometries as well as present 
applications for complex targets, 

• demonstrate stabilization of EFIE-IGA solvers using the proposed basis sets together with 
frequency scaling, 

• demonstrate construction of the Galderon preconditioner using the proposed basis sets (with¬ 
out the need for auxiliary barycentric meshes), 

• and, present scattering data from multiscale obstacles. 

The paper is organized as follows: In Section [^outlines the problem, whereas Section|^presents 
a brief summary of extant literature on subdivision surfaces. Section details the crux of this 
paper: (i) defines basis functions on the subdivision surface, (ii) develops methods to address low 
frequency breakdown and multiplicative Calderon preconditioners for the EFIE. Section [^presents 
several numerical results that verify the accuracy and efficacy of this approach. Finally, Section]^ 
summarizes the contribution of the work as well as outlines directions for future research. 


2. Integral equations for Electromagnetic Scattering 

The model problem for analysis will be a perfect electrically conducting object that occupies 
a volume whose surface is denoted by dfl. It is assumed that this surface is equipped with an 
outward pointing normal denoted by h(r). The region external to this volume (IR^\r2) is occupied 
by free space. It is assumed that a plane wave characterized by {E*(r), H*(r)} is incident upon this 
object. The scattered field for r G can be obtained using equivalence theorems as follows: 

h(r) X E®(r) = To J(r) 
h(r) X H®(r) = /C o J(r) 


where 


To J(r) 


-JWMoh(r) X 




/Co J(r) 


h(r) X V X 



dr'g{r, r')J(r) 


(lb) 


where g{r,r') = exp [— jKjr — r'|]/(47r|r — r'|), k is the wave number in free space, J(r) is the 
equivalent current that is induced on surface, and I is the idempotent. In the above expression 
and in everything that follows, we have implicitly assumed (and suppressed) an exp [jujt] time 
dependence. Using the above equations, one may prescribe the requisite integral equations as 


EFIE: n{r) x h(r) x (E*(r) -|- E®(r)) =0 Vr G 
MFIE: h(r) x (H*(r) -f H®(r)) = 0 Vr G dQ- 


(2a) 

(2b) 


that are known as the electric/magnetic field integral equations (EFIE/MFIE). In (2b), de¬ 
notes a surface that is conformal to but just inside dVl. These equations do not have unique solutions 
at the so-called irregular frequencies, but may be combined to yield the combined field integral 
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equation m that gives unique solution at all frequencies. It should be noted that while these are 
the most popular formulations that are used in practice, they are not the only ones. Other formu¬ 
lations such as the combined source integral equation [28] , augmented EFIE, augmented MFIE [29] 
and the charge current integral equation |30j exist and have seen recent development. However, the 
bulk of this paper’s focus will be the EFIE and its discretization. The choice is largely motivated by 
the numerous challenges that exists in solving these equations, both at the mid and low frequency 
regimes. To this end, in what follows, we will develop basis functions that are well formulated and 
present modifications to the formulation to account for both low frequency behavior and to impose 
Calderon preconditioners. 

Methods for solving these equations follow the standard prescriptions: (i) represent the surface 
of the scatterer using some tessellation, (ii) define basis function on this approximation to the 
geometry, (iii) demonstrate desirable properties of these basis sets, and (iv) validate solutions 
to integral equations solved using this procedure. As elucidated earlier, as opposed to classical 
tessellation, we use subdivision surfaces for geometric representation. This is elucidated next. 

3. Box Splines and Subdivision Surface 

In what follows, we provide a brief overview of shape description as effected by subdivision 
surfaces. This section is provided purely for completeness and omits details that can be found 
in several references, e.g., jSIl EH ESI- The traditional workhorse of geometric description is 
NURBS. As a result, current isogeometric methods have been largely developed using NURBS 
as basis functions. However, a singular feature stands out that poses to be a bottleneck-while 
NURBS descriptions are in the interior of a patch, the continuity may be C® or even worse 
(not watertight) at the boundary between patches or at intersection curves. As a result, one has 
to define physical basis functions that are div- or curl-conforming for vector problems. This stands 
in contrast to subdivision surfaces that are constructed as limit surfaces obtained by subdivision 
processes. This leads to meshes that are almost everywhere except at isolated points where 
they are C^. 

As triangular tessellation is ubiquitous in computational electromagnetics, the choice of the 
presentation below is based on the Loop subdivision scheme for closed triangular surface meshes. 
Let Mq denote the primal base mesh that leads to the limit surface. As in any tessellation, Mq 
comprises a set of vertices V and a connectivity map. The 1-ring (union of incident triangles) of a 
vertex produced by the connectivity map can be characterized by the valence of the vertex; specifi¬ 
cally, (i) a valence-6 node/vertex is deemed regular and (ii) any other vertex is called extraordinary 
or irregular. A triangle is regular if all its vertices are regular, and irregular otherwise. Fig. 
illustrates an irregular triangle and one subdivision around it. 



Figure 1: Irregular vertex and triangle (a) An irregular triangle (vertex 1 is valence-7), 
(b)subdivision once [33] 
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Figure 2: Local averaging during refinement. 


The Loop subdivision scheme defines a smooth surface approximating positions assigned to each 
vertex in the base mesh through a procedure (Fig. [^. Starting from the base mesh, an infinite 
sequence of meshes is determined through repeatedly refining the fc-th mesh into the (fc+l)-th mesh 
by dividing each triangle into four subtriangles by inserting one vertex per edge, and adjusting 
vertex locations. The location of the existing vertex on the refined mesh is an affine combination 
of its one-ring vertices with weight = (5/8 — (3/8 -I- cos(27r/n)/4)^)/n, and itself with weight 
1 — na„, where n is the valence. The location of the newly inserted vertex is an affine combination 
of the end vertices of the edge, each with weight 3/8, and the other two vertices of the two incident 
triangles, each with weight 1/8. The limit of the sequence of piecewise flat surface is the Loop 
subdivision surface with smoothness almost everywhere, except at the extraordinary vertices, 
where it has smoothness. Note, the newly inserted vertices will always be regular (valence-6). 
Thus, we can assume that extraordinary vertices are at least two edges away from each other, 
otherwise the base mesh is replaced by the mesh obtained after proceeding through the subdivision 
process once. 

The same procedure through repeated refinement with local averaging can be used to define 
smooth functions on the Loop subdivision surface. In the context of isogeometric analysis, it is 
perhaps useful to view this process in terms of effective basis functions. Consider the a limit 
surface S(u,v) = that is defined by coordinates of the vertices Vi and ^i{u,v) is an 

effective basis function that is associated with quantities associated with vertex 1/, and {u,v) are 
the pairwise coordinates on a parameterization chart. For instance, ^i(u,v) can be constructed by 
the subdivision procedure starting with 1 assigned to Vi and 0 assigned to all other vertices. 

Such a basis function ^i{u, v) has a compact support, as it is 0 outside C^’s 2-ring (union of the 
1-rings of the 1-ring vertices). Piecewise-linear parameterizations of the base mesh can be used as 
a parameterization of the subdivision surface. In a parameterization chart (m,u), it can be shown 
that for any vertex Vi, ^i{u, v) can be expressed in closed form as the 3-direction quartic box spline 
Ni{u,v) if {u,v) is located within a regular triangle. To illustrate this, for triangles incident on 
Vi (triangles with label 3 in Fig. [^, e.g., triangle denoting the linear basis function (hat 

function) of Vi by (pi(u, v) and setting r = ipi, s = ipj and t = ipk (note that r -|- s -I- t = 1), 

Ni{u, v) = 6r^ -I- 24r^t -|- 24r^/^ -|- 8rt^ +1"^ + 24r^s -|- 60r^s/ -I- 36rst^ 

+6st^ + 24r^s^ -I- 36rs^t + \2s^t^ + 8rs^ -I- Qs^t + s'*. 

For triangles with 2 vertices in the 1-ring of Vi (triangles with label 2 in Fig. e.g., triangle T^ik, 
setting r = ipj, s = (pi and t = fk, 

Ni{u, v) =r''^ + 6r^t + -|- 6rt** + f + 2r^s + 6r^st + 6rs/^ -I- 2st^. 

For triangles with one vertex in the 1-ring of Vi (triangles with label 1 in Fig. [^, e.g., triangle 
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Tjmi, setting r = ipj and t = (pi, 


Ni{u, v) = + 2r^t. 



Figure 3: 2-ring of a valence-/c vertex (A: = 9 in the example). 


For an irregular triangle Tijk with Vi being extraordinary, the evaluation of for Vi in the 
one-ring of Vi, Vj or 14 is modified. If pi(u,v) = 1, the evaluation is right at the extraordinary 
vertex with valence n, where it can be shown that = l/(g;^—h n) for V) in t^i’s 1-ring, 

and ^i{u,v) = 1 — n/(g|—h n). Otherwise, when 2“^ < ipi{u,v) < the point belongs to 

Oj in Fig. 1^ if (pj{u,v) > 2“^', to Og if ipk{u,v) > 2“*, or to otherwise. It means after k-th 
subdivision, the basis functions (^i, and ^k) can be evaluated through linear combinations of box 
splines with values on the relevant vertices after k-th subdivision, since is a regular triangle 
after fc-th snbdivision. The relevant values for each can be stored in a vector Ek^m containing 
12 scalars, one for each vertex whose box spline basis function has the triangle in its support, which 
can be calculated as 

Ek,m = PmAA^ Eq, 

where Eq is the initial assignment of the n -|- 6 vertices in Fig. |la[ which is 1 for V) when evaluating 
, and 0 for any other vertex; A is the subdivision operator matrix, an (n -|- 6) x (n -I- 6)-matrix 
denoting how the values influencing the irregular triangle in the next level are calculated from the 
values at the current level (Fig. lb shows the vertex at the next level); A is an extended version 
of A including how the vertices n -\- 1 through n -I- 12 of next level is computed form the current 
level; Pm is the picking matrix selecting the relevant 12 values for within the n -I- 12 values. 
Note that, in practical implementation, eigenanalysis of A is performed so the evaluation involves 
taking powers of the eigenvalues as an efficient way of handling the calculation of A^. Furthermore, 
the subdominant eigenvectors provide a means for direct evaluation of derivatives even in irregular 
triangles [33]. 

Several properties of these basis functions are worth noting as they are critical to represent 
both the geometry as well as the physics defined on the geometry. Note, some of these properties 
are available for a NURBS description as well in the interior of the element, but not necessarily 
across domains. These include the following: 


1. compact support, 

2. non-negativity, 

3. convexity preserving, 

4. and continuity across patch boundaries. 

The above properties make this subdivision basis functions a good candidate for modeling physics 
on both regular and irregular (with the help of subdivision scheme) triangular meshes. 
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Figure 4: The parameter domain within an irregular triangle is recursively partitioned into an 
infinite number of subdomains [33]. 


4. Current representation and field solvers 

4-.1. Isogeometric basis sets 

Section developed a framework wherein one dealt with representation of the limit geometry 
surface. Using the same representation as presented earlier, we introduce the concept of scalar 
functions that are defined on the limit surface. As with geometry, the “limit” function can be 
represented in terms of 1-ring control weights of a regular triangle and the box splines as 

12 

/(r) = f{r{u, v)) = ^ N,{u, v)wi\ for r e T (3) 

i=l 

In the above, we have assumed that the mapping r(u,v) exists, where (u,v) are the barycentric 
coordinates of a triangle. In what follows, this assumption will be implicit, and only r will be used. 
It can be readily verified that if only one vertex has a weight of unity and other zeros, one would 
immediately get a smooth (up to 2nd order continuity globally) scalar function. As a result, one 
can associate an effective basis function with every vertex such that the limit function 

Nv 

/(r) = (4) 

n—1 

where Ny is the number of vertex and ^n(r) is an effective basis function that describes the 
influence of the scalar quantities associated with a vertex. Fig. gives an example of the scalar 
basis function used for formulating isogeometric analysis on top of subdivision-described surface. 
The basis function is associated with an irregular vertex of valence 8. To demonstrate 2nd order 
continuity, the surface Laplacian of the scalar function is plotted in Fig. |5b| The behavior of 
the subdivision basis and its derivatives follow: ^n(r) ~ 0{1), Vs^Ti(r) « 0{l/h) and Vg^„(r) « 
0{l/h^) where h is an approximate dimension of the patch. 

Properties of representation. Two salient properties of the scalar basis used for this definition are: 
(i) as they rely on approximating subdivision they are (7^ almost everywhere; and (ii) if U = 
where is the domain associated with ^n(r), then the function ^„(r) vanishes on In other 

words, the basis functions ^n(r) G C§ almost everywhere. 

Next, to use this basis set to represent the current, we note that current on any surface can be 
represented using a Helmholtz decomposition as follows: 

J(r) = Vs(()(r) + V X (n'0(r)) -I- n7(r) (5) 
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(a) (b) 

Figure 5: Scalar basis function associated with a vertex with valence of 8. (a) Basis function, (b) 
Surface Laplacian of the basis function 


where 'n7(r) is the harmonic field and has a trivial solution (ra7(r) = 0) for simply-connected struc¬ 
tures, and (j){r) can be associated with the charge space. Following the subdivision representation 
in Q, we will assume that ^(r) and '(/j(r) represent approximations to ^(r) and ?/'(r), respectively, 
and they can be represented in a manner similar to Q; viz., 

^(r) « ^(r) = ^ai,„Cn(r) 

( 6 ) 

V'(r) « ^{r) = ^a2,„Cn(r) 


Using Q, it is now possible to define the approximations to the current as 

J(r) « J(r) = ^ (r) -b a2,„J^(r)] 

n 

Ji = \/sUr) 

Jl = hx VsC„(r) 


(7) 


The physical interpretation of this representation is akin to a standard subdivision representation 
of a limit surface; we represent the the “limit” current via the “limit” scalar functions. We note, 
that the basis for the currents (and the auxiliary potentials) are represented via operations on the 
subdivision basis which can be effected numerically rather trivially. In this work, the standard 
definition of the inner product (X, Y) = drX ■ Y is used. 


Properties of the basis function. Several properties of the basis functions make this definition 
appealing. These are as follows: 

Continuity Due to the reliance on functions 'Cn(r) that are Cg almost everywhere, the resulting basis 
functions have continuity almost everywhere and continuity at isolated points. This stands 
in stark contrast with classical Rao-Wilton-Glisson (RWG) or their higher order counterparts 
that are div-conforming but not at the boundary. 

Orthogonality The inner product of (J)j,J^) = 0. The proof for this assertion can be trivially 
obtained using Green’s theorems together with properties of the surface curl and the fact that 
Cn(r) = 0 for (r) s 9D. Specifically, 


f drJlir) ■ 3l{r) = (f drf„{r)u ■ 3l{r) - f dr^„(r)V^ • J^(r) 
jQ,n JdQn J 


= 0 


( 8 ) 


where u is the normal of tangent to the surface. 

Mapping The function J^(r) maps only to Vs^(r) and J^(r) maps only to h x Ws'fp{r). The proof 
can be obtained using Green’s theorems, properties of the surface curl/divergence, and the fact 
that i{u,v) is a partition of unity. Specifically, 



f dr^„(r)u • J(r) 
Jdn„ 



dr^„(r)Vs • J(r) 



dr^„(r)V2(^(r) 


(9) 


The key take home message is that J^(r) completely represents the charge space and the rest of 
the current is represented by J^(r). 

Charge neutrality The basis functions maintain charge neutrality. This is proved using 


drVs • J(r) = V ai,„ / drVs ■ VsC„(r) 

! „ Jn„ 

= y'ai,n'P drii ■ V 

„ ddn„ 


( 10 ) 


= 0 as Vs^n(r) = 0 Vr s dflr, 


Compactness The definition of basis function is local. 

Refinement The basis functions are subdivision based. As a result, they inherit properties of subdi¬ 
vision representation including adaptivity. 

These properties ensure a complete representation of the currents on the surface of an simply 
connected object, and forms a rigorous Helmholtz decomposition of currents on the surface. Again, 
as opposed to classical basis functions defined on tessellation, Helmholtz decomposition is inherent 
in the definition of basis functions where as the common approach is to design “quasi”-Helmholtz 
decomposition using weighted sum of basis functions. It should be noted that, in the above rep¬ 
resentation, the potential functions are the actual unknown functions and not the currents. As a 
result, extra steps have to be taken to ensure their uniqueness. These are elucidated below when 
using these basis functions within a Galerkin framework. 

Using these definition of basis functions, each vertex is associated with two degrees of freedom. 
Fig. [6] shows the behavior of the two basis functions associated with a node of valence 8. A point 
to note is that the number of solenoidal and non-solenoidal basis functions are equal to each other. 
As we will see, this helps creating Galderon preconditioners. 


4-S- Field Solvers 

Given the prescription of basis functions, the discretized version of (2a I can be obtained using 
Galerkin testing. The resulting matrix equation can be written as 


■ Zii 



(11a) 


where 


^nm — j^d'O / 

drJ(^(r) 

• / dr'g{r,r')3’^{r') 


1 

/ drVs ■ 

Jnm 

Jn)!") / '^r' 5 (i’,r')Vs ■ J)„(r') 

(11b) 

WEo J 

n„ 

Jam 


= 

-^n 


= [ drJ^(r).E*(r) 

(11c) 


In the above equations, Sij denotes a Kronecker’s delta. Interesting features of the above expression 
are apparent; (i) is the only term that has the charge contribution. As will be shown later. 
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(a) (b) 

Figure 6: Basis function associated with a vertex with valence of 8. (a) nonsolenoidal type, (b) 
solenoidal type 


this “decoupling” is an essential component for the construction of low frequency stable solvers, 
and (ii) since the system of equations are constructed using conditions on currents that rely on 
derivatives of the potential 0(r) and ^/’(r), it follows that the potentials are determined upto a 
constant, and the number of degrees of freedom is less by one for each potential. While one can 
impose this via different means, we have chosen to essentially constrain the system by choosing 
ai^N and 02 ,at to be zero. This implies a trivial change to the system of equations ( [Tl] ). The 
evaluation of the inner products in the above equations is effected via higher order quadrature 
and Duffy integration rules. While the approach presented thus far is valid for all frequencies, it 
suffers from low frequency breakdown. In what follows, we demonstrate that (111 can be modified 
trivially to alleviate low-frequency breakdown. Likewise, Calderon preconditioners can be easily 
implemented using this space of basis functions. 

But before we proceed, interesting insight may be gained by examining the inner products 
that arise using the Galerkin procedure. Matrices arise from testing the electric field, i.e., 
corresponds to measuring the radiated electric field due to Em(r) = E {J^} by J^. If follows that 



f c?rJ^(r) • E^(r) = /" drVs^„(r) • E^(r) 

= f dr^ri{r)u ■'Em{r) - dr^„(r)Vs • Em(r) 

Jan„ 

= - dr^„(r)Vs • E„(r) 

Jfi„ 

drJ^(r) • Em(r) = [ drh x VsCn(r) • E™(r) 

= [ drn X VC„(r) • E„(r) 

= [ drn- (V^„(r) x E^(r)) 

= [ drn • V X (^„(r)E„(r)) - /" dr^„(r)n • V x E^(r) 

= f dr^„(r)Em(r) • f-b jw/io / dr^„(r)n • Hm(r) 

= ju}fj.o / dr^„(r)n • Hm(r) 

Jn„ 


(12a) 


(12b) 
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In the above equation, i is a unit vector tangential to the boundary dilm and H„(r) is the magnetic 
field due to J^. From the above equations, it is apparent that (12aI tests the tangential component 
of the electric field. However, (12b) yields equations that test in normal component of the magnetic 
field. So the two basis functions used in the analysis impose tangential continuity of the electric 
field as well as the normal component of the magnetic field. As a result, using a basis that satisfies 
the Helmholtz decomposition (and relies on scalar function that are Cq), results in equations that 
naturally fit into the framework of the Current-Charge Integral Equations (CCIE) [3D]. It should 
also be noted that while eqns. (12a) and (12b) were specified for the scattered field, they are 
equally valid for the incident field. Indeed, using the final expression in to evaluate the 

integrals is more accurate and requires fewer quadrature points. 


4-3. Low-frequency stable EFIE 

It is well known that the EFIE suffers from low-frequency breakdown. The rationale for this 
breakdown is readily apparent by examining the components of the impedance matrix in (Hill- 
Assuming that the average largest linear dimension of the support of the patch is h. It follows that 
the entries corresponding to 

Z't = 0{^h^) + SnSkiOil/n) (13) 


These results follow from the fact that the functions ^(r) = 0{1). The above scaling indicates 
that a portion of the elements associated with source and test basis functions that are associated 
with irrotational functions J)i(r) scale as 0{1/k), whereas all others scale as 0{Kh'^). This implies 
that as K —^ 0, the portion of that corresponds to the charge contribution dwarfs the rest. 
This situation is similar to those encountered for by classical Nedelec elements. Here, one takes 
recourse to loop-star/loop-tree decompositions |3H|33 |3D] that effect an approximate Helmholtz 
decomposition of the currents. As has been shown, the resulting decompositions contain a portion 
that is exactly divergence free and one that is approximately curl free. Whereas the support of 
the divergence-free portion is local, the same is not true of curl-free portion. This is in contrast to 
the method presented here wherein both components have local support. Using these basis results 
in the charge being modeled correctly, and a matrix whose scaling looks like (13). Rescaling these 
equations has been shown to render the solution stable. Following a similar procedure, it can be 
shown that rescaling both the matrix elements, the coefficients and the right hand side results in 


where 


= jujl3iktio 

[ drJlir) ■ 

/ d.r'g{r,r')3^{r') 




jSiiSki 1 

^ drV. • J), (r) / dr'g(r,r')V«-Ji,(r') 

£o J 





if 1 = l,k = 1 

Pik 


if ; = 2, fc = 2 


1 1 

otherwise 



/ drJJ((r) . E^r) 


(14) 


(15) 

(16) 


As prescribed, these equations achieve the desired stability. As can be trivially shown, these equa¬ 
tions decouple as w —> 0. This is akin to similar prescriptions for classical loop-star/tree algorithms 
|34l|35l|37]. However, a salient feature of the IGA-basis is that constructing a well behaved system 
is tantamount to using a diagonal preconditioning sans constructing the complementary system 
that is required for a loop-tree/star (Hodge) decomposition. 

It should be noted that Helmholtz decomposition is not the only way to solve low frequency 
breakdown, and other techniques 11113 GB Emi nni Gi] exist. The isogeometric basis functions pre¬ 
sented herein could help construct those methods. 


4.4- Calderon preconditioner 

It is well known that the standard EFIE presented in the earlier sections is a first kind equation, 
and as with all first kind equations, can be ill-conditioned especially when the spatial scales in the 
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problem are widely separated. As has been shown by several others (see Refs. mum mi, and 
references therein), the rationale for the ill-conditioning is the spectral separation of the eigenvalues 
of the operator with increase in discretization density, with a set that clusters around zero and 
others that cluster around infinity. Given this separation, the resulting matrix systems become 
rapidly ill-conditioned. The remedy to this problem exploits the Calderon identities wherein the 
EFIE operator (T {•}) preconditions itself resulting in a second kind integral operator whose eigen¬ 
values accumulate around —1/4 [32113 El]- The challenge in using these identities was the lack 
of well behaved Gram matrices that link the domain and range of the T operator|46l [47]. For 
Thomas-Raviart/Rao-Wilton-Glisson basis sets and (7° geometries, basis functions developed by 
Buffa and Christiansen [35] have been exploited in a sequence of papers to thoroughly understand 
and solve this problem. In what follows, we show that basis functions defined herein result in a well 
conditioned Gram matrix, thereby permitting a natural discretization of the Calderon operator. 
To begin, we define operators 

roJ(r) =raOj(r)+r^oJ(r) (17a) 

where 


77 o J(r) =-jw^o^(r) X / dr'g(r, r') • J(r') 

Jan 

7/, o J(r) = —J —— h(r) X f dr'Wg{r,r') ■ 
^^0 Jdfl 


(17b) 


The Calderon projector is defined as 7"^ o J(r) = To T o J(r). It has been shown that o 
J(r) = (—1/4 -I- A^) o J(r), where K? o J(r) = K. o jC o J(r) is compact, and the eigenval¬ 
ues are clustering around —1/4. Alternatively, this operator may also be written as o J = 
(77 077+77077 + 77 , 077+770 77 ) o J. Since direct discretization of is impossible, typical 
implementation follows the multiplier approach; i.e., define intermediate mapping from the range 
space of the T operator to the domain of the T. This is typically effected via a Gram matrix and 
has been extensively explored. Given the usage of Helmholtz decomposition Q and the orthogo¬ 
nality ([^ between the two components, it can be shown that 77 o 77 o J = 0. Enforcing the latter 
condition has been difficult to accomplish in traditional discretization and representation schemes. 

The principal advantage of using Q is that it enforces an exact global Helmholtz decomposition 
in terms of a set of local scalar functions that are Cq almost everywhere. Recall the following: 
(J7, Jm) = where is the result of evaluating the inner product over the support of 

basis functions H 0^ and is defined in (19). As a result, the Gram matrix Q is block diagonal 
with entries = InmSik- Given that, for a pair of spaces, (J^, J^), it can be shown that the 

following sequence is satisfied: (J^,J^) ^ (J^,J^) ^ (J^ - t 2 \ T, 

J7P) 


and (J^, J^) TT- (J^, J^). It follows from the above that 
form, the above sequence can be rewritten as 


TQ-^T = 


>J^)- 
(J7J2' 


Specifically, (J^,J^) —(J^,0) 
(0,0)- In matrix 


r 

r2i -1 

• a 

■ 

0 

-1 

r j-ii 

' a 

r2i -1 

• a 


q~22 

'a 

0 

^22 


n-21 , n-21 

1 f a ^ f(j) 

q-22 

a 


(18) 


Finally, a critical component in discretizing the T^ operator is the Gram matrix. The benefits of 
using local Helmholtz decomposition becomes apparent in its construction with the off-diagonal 
blocks being zero and the element of the diagonal blocks being 


lnni= / drVsin{Y) ■ Vs^m{Y) (19) 

jQrinQm. 

that are a variational form of the Laplace-Beltrami operator on a scalar function. The resulting 
system is positive definite leading to well behaved inverse of the Gram system. 


5. Numerical Examples 

This section presents several numerical examples to demonstrate the efficacy of the proposed 
approach to electromagnetic analysis. In order to do so, we shall present data that demonstrates 
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Table 1: Rel. L 2 error in the currents density 


subdivision times 

0 

1 

2 

average element size 

0.1511 

0.0755 

0.0377 

Rel. L 2 error 

0.0257 

0.0067 

0.0030 


the following; (i) accuracy of the approach to when compared against analytical data; (ii) examples 
illustrating low frequency stability, (iii) analysis of multiscale structures to illustrate the efficiency 
of the Calderon preconditioner, and (iv) application to complex targets,. Unless stated otherwise, 
the data presented compares in radar scattering cross-section computed using the IGA-MoM solver 
and those computed using either analytical methods (if available) or a validated method of moments 
code that is based on RWG basis functions. 

5.1. Accuracy of IGA-MoM 

To validate and demonstrate the accuracy of the proposed approach, we consider scattering 
from a sphere discretized at multiple resolutions. To this end, consider a sphere with radius of lA 
that is modeled using an initial control mesh comprising 642 vertices. Starting from this control 
mesh and performing one and two Loop subdivision results in two other meshes with 2562 and 
10242 vertices, respectively. This results in three meshes with 1280, 5120, and 20480 faces. Note, 
the limit surface for all three mesh densities is identical. As a result, what changes is the support 
of the basis functions, and therefore, approximations to the discrete operators. In all the data 
presented below, a field propagating in fc = —z and polarized along x axis is incident on the 
sphere. The number of degrees of freedom (DoF) for IGA-MoM is twice the number of vertices; 
consequently, the three tests have 1284, 5124, and 20484 DOFs, respectively. Table presents 
errors in currents on the surface of the sphere (between those obtained using IGA-MoM and Mie 
series) for these three different discretizations. Likewise, Figs. and present pointwise errors in 
absolute values of the real and imaginary parts of the current for the case of two subdivisions. As 
is evident, the accuracy in the current both in the L 2 and Li norms is excellent. 



(a) (b) 

Figure 7: Magnitude of surface currents density on the sphere: (a) real part and (b) imaginary 
part 


5.2. Scattering from Structures at regular frequency 

In what follows, we analyze the performance of the proposed approach and compare these with 
a conventional MoM EFIE solution technique that relies on RWG basis functions. Unless specified 
otherwise, the metric for comparison is RGS data in the (/) = 0 plane. Eurthermore, the geometry 
of all scatterers that are analyzed in this paper is smooth; extension to non-smooth structures is 
non-trivial but realizable from a subdivision perspective [49] . IGA solvers that include in open and 
sharp surfaces will form the basis of a future work. 

The first example is a truncated cone, with the height 3.0A and the radii of the top and bottom 
circular cross-sections being 0.2A and lA, respectively. A plane wave that is polarized along x and 
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Figure 8: Radar cross section of the sphere {(j) = 0 cut) 




Pseudocolor 
Van abs. 
r 2.954e-05 
-2.215e-05 


Max: 

Min: 


Max: 

Min: 


Pseudocolor 
Var: val 

4.482e-05 
3.362e-05 




(a) 


(b) 


Figure 9: Pointwise relative error (a) real part, (b) imaginary part with 20484 DoFs 


propagating in the —z direction is incident on the object. The object is represented using 12378 
nodes; this number is increased at the top and bottom surfaces so as to maintain a sufficient degree 
of sharpness. The number of DoFs for IGA-MoM and conventional MoM are 24756 and 37128, 
respectively. The surface current density obtained using IGA-MoM is depicted in Fig. It is 
evident from these plots that the currents are smooth (without unphysical aberrations anywhere, 
especially near the crease). Fig. 11 compares the radar cross section between the proposed method 
and the conventional method of moments, and it is apparent that the agreement between the two 
is excellent. 

The second example is a structure composed of a block (3.33A x 10.OA x 1.33A) and five cylinders 
(radius is 0.5A and the height is 0.67A) uniformly distributed on the top surface of the block. An 
electromagnetic held that is propagating along —y and polarized along the z incident on this object. 
The number of DoF for the IGA-MoM and conventional MoM are 29028 and 43536, respectively. 
Figur^l^ presents the real and imaginary parts of the current distribution on the surface of the 
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Pseudocolor 



Pseudocolor 
Var: Imoginary 




(a) (b) 

Figure 10: Magnitude of surface current density on the truncated cone: (a) real part and (b) 
imaginary part 



Scattering angle/deg 


Figure 11: Radar cross section of the truncated cone (0 = 0 cut) 


object, and again, it is evident that the results are smooth without artifacts. Further, excellent 
agreement can be seen between IGA-MoM and conventional MoM in the RCS data. 

As a final example in this set, we analyze scattering from a model airplane whose dimensions are 
5.49A X 5.48A x 1.52A. The plane wave incident on the object is propagating along —y and polarized 
along z. As before, data is obtained using both IGA-MoM and conventional MoM using 23988 and 
35976 DoFs, respectively. The current distribution using IGA-MoM is depicted in Fig fM} Again, 
it is evident that the current is well captured as is the scattering cross-section (see Fig. |15|). 

5.3. Conditioning of Isogeometric System at Low-frequency 

Next, we discuss the stability of the resulting system (and its preconditioned form) at low 
frequency using IGA basis and the aforementioned wavenumber scaling. The metrics used for the 
discussion are both the condition number of the resulting system and the number of iterations 
necessary to achieve a desired error tolerance when solving this system using an iterative solver. 
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Pseudocolor 


(a) 


Pseudocolor 
Var: Imoginary 



(b) 


Figure 12: Magnitude of surface current density on the composite structure: (a) real part and (b) 
imaginary part 



Scattering angle/deg 


Figure 13: Radar cross section of the composite structure {(j) = 90 cut) 


In all the examples presented, we use a generalized minimal residual (GMRES) iterative solver 
and the object being analyzed is a sphere with radius Im that is discretized such that the average 
distance between vertices is around 0.15m. 

Fig. depicts the condition number of IGA-MoM systems as a function of frequency from 
IHz to lOOMHz. The three curves are, respectively, for the original IGA system, the diagonal 
preconditioned system, the Galderon preconditioned system. From the plots, it is evident that 
the condition number is almost constant over the whole band for all the systems. Both diagonal 
preconditioner and application of the Galderon precondition offer significant improvements over the 
original IGA-MoM. Note, the behavior of the original IGA-MoM is very much unlike conventional 
MoM, thanks to the properties of the basis function used. 

Next, we examine the number of GMRES iterations required to converge to an error tolerance 
of 1.0 X 10“® and 1.0 x 10“^° at various frequencies. These results are depicted in Eigs. 17 and 
Fig. 18 respectively, both sampled at IHz, lOOHz, lO^^Hz, 10®Hz and 10®Hz. Since the iterative 
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(a) (b) 

Figure 14: Magnitude of surface current density on the air plane model: (a) real part and (b) 
imaginary part 



Scattering angle/deg 


Figure 15: Radar cross section of the air plane model {(p = 90 cut) 


numbers for the original IGA system are much higher than that of the diagonal preconditioned and 
Calderon preconditioned system, only the iteration numbers for the latter two systems are given. 
As seen in the two plots, the number of iterations required for both systems are relatively constant 
across a range of frequencies (from IHz to lOOMHz). 

These two tests serve to illustrate some of the salient features of the proposed approach, (i) 
Low frequency stability is achieved simply by diagonal preconditioning. This is in contrast with 
conventional approaches where it is explicitly necessary to define auxiliary unknowns albeit at linear 
cost, (ii) Imposition of the Calderon preconditioner involves an inversion of the Gramm matrix 
as an addition operation. The original matrix system is retained which makes integration with 
fast methods trivial. This is different from existing methods that require Buffa-Christiansen basis 
functions to be defined on barycentric meshes. Finally, as is evident from Fig. the results 
obtained from the Calderon preconditioned solver at IHz agrees very well with analytical data 
obtained. 
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Figure 16: Condition number at different frequencies 
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Figure 17: Number of GMRES iterations to converge to a tolerance of 1.0 x 10 ® 
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Figure 18: Number of GMRES iterations to converge to a tolerance of 1.0 x 10 



Figure 19: Comparison of RCS data obtained at IHz with analytical solutions 
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5.4- Examples with Multi-scale mesh 

In all the examples analyzed thus far, the initial control mesh is such that the resulting tessela- 
tion is almost uniform in that the ratio of the maximum edge length to the minimum edge length 
is 0(1). However, a more intellectually interesting and practically applicable problem is when this 
ratio is significantly higher. It is challenging to design stable methods for these problems due to two 
effects that act in concert with each other; wavenumber and element size scaling. In this example, 
scattering from a sphere with locally refined meshes is simulated, and verified by comparing with 
the one with almost uniform initial control mesh. Fig. [^demonstrates the multiscale control mesh 
for the unit sphere. The radius of the sphere is lA, and the incident plane wave is polarized along 
the X direction and propagating in the —z direction. The resulting system of equations is solved 
using GMRES. Figure compares the convergence history for both the Calderon preconditioned 
system as well as a diagonal preconditioned system. The efficacy of the Calderon preconditioned 
system is evident by the rapid convergence. 



(a) (b) 


Figure 20: Multiscale control mesh: (a) the whole mesh and (b) locally refined region 


6. Summary 

In this paper, we have developed a novel isogeometric analysis technique for solving the electric 
field integral equation that is encountered in electromagnetic field analysis. The fundamental 
thesis of isogeometric methods is that they inherit the rather signihcant advantages of modern 
CAD tools-smoothness of geometry representation, morphing, dynamic meshes (if necessary), etc- 
by using the same basis sets for both geometry representation as well as representation of physics 
on this geometry. The choice of the electric field integral equation is predicated upon the fact 
that it is one of the more challenging equations to solve. The approach presented here leverages 
existing geometric construction techniques that use subdivision to define basis functions that result 
in well behaved integral operators. Thanks to these operators, it is possible to trivially modify 
these to impose Calderon preconditioners, and construct systems that are low frequency stable 
and can handle multiscale geometric features. A number of results demonstrate the convergence 
and accuracy of the technique, applicability to computation of scattering at both regular and 
low frequencies, as well as to structures with multiscale features. It is very simple to apply the 
proposed basis set to other types of integral equations, such as magnetic field integral equation and 
combined held integral equation. However, as with the introduction of any new technique several 
open problems remain: these include addition of features to handle edges and open domains, 
integration with fast solvers, development of techniques for multiply connected objects, using this 
framework to include Debye sources, etc. Several of these topics are active areas of research within 
the group and will be presented in different forums soon. 
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Figure 21: Convergence of the iterative solver 
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